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A MODIFIED METHOD OF INTEGRAL RELATIONS APPROACH TO THE 
BLUNT-BODY EQUILIBRIUM AIR FLOW FIELD, INCLUDING 
COMPARISONS WITH INVERSE SOLUTIONS 

By L. Bernard Garrett, John T. Suttles 
Langley Research Center 

and 

John N. Perkins 
North Carolina State University 

SUMMARY 

A modified approach to the first-order approximation of the method of integral 
relations is presented for the numerical calculation of the inviscid, adiabatic flow field 
around a blunt-nose body traveling at hypersonic speeds. Solutions have been obtained 
and computed flow-field data are presented for the adiabatic steady flow of a perfect gas 
and air in chemical equilibrium. The results obtained by the present method of solution 
are compared with results obtained by using inverse methods, property-derivative inte- 
gral relations methods, and experiment. The results indicate that the present method 
provides an accurate, versatile solution of the blunt-body flow field, and the relative sim- 
plicity of the method should allow extension to coupled radiating flow-field analyses. 

INTRODUCTION 

Design for interplanetary spacecraft depends on a realistic definition of the entire 
flow field surrounding the body. The flow-field solutions should provide the necessary 
inputs required for the calculation of aerodynamic loads, shear stresses, radio attenua- 
tion regions, and convective and radiative heating. At hyperbolic speeds a predominant 
consideration is the radiation fluxes which are strongly dependent upon the distribution 
of the fluid properties and the chemical species between the shock and the body. At these 
velocities the energy losses due to radiation can have a significant effect on the entire 
flow field. 

The purpose of this report is to develop an approximate numerical method which 
has promise for application in coupled radiating flow-field analyses and to compare 



adiabatic-shock-layer results with other existing numerical solutions and with experi- 
mental data. 

A method of integral-relations approach, which was first applied by Belotserkovskii 
(ref. 1) and Traugott (ref. 2) to the flow of a perfect gas over a blunt body, is employed 
in this analysis. A one-strip approximation is used. The solution is obtained by inte- 
grating a set of governing differential equations which are cast in terms of the "fluxes'' 
(mass, momentum, etc.) or "complexes" which was originally recommended by Lun'kin 
et al. (ref. 3) and Kao (ref. 4) and is employed extensively by Belotserkovskii (ref. 5) 
near the sonic region of the flow field. This approach is in contrast to prior approaches 
(see ref. 6, for example) in which the governing equations are written in terms of deriva- 
tives of properties (velocity, pressure, etc.). The flow equations are developed without 
imposing the constancy of the entropy along streamlines; thus, the method can be extended 
to coupled radiation flow-field analyses. The approach is made more flexible in that the 
thermodynamic characterization of the fluid does not enter directly into the governing 
differential equations, but enters the governing system of equations only through auxiliary 
equations of state. 

Flow-field results are presented for the adiabatic steady flow of a perfect gas and 
air in chemical equilibrium. The primary goal of this research is to examine the equi- 
librium air results obtained from the relatively simple one-strip integral approach to 
determine whether the solution is sufficiently accurate for use in adiabatic radiation com- 
putations such as were done by Kuby, et al. (ref. 7). 

The perfect gas model is employed primarily to study the basic flow-field approach, 
namely, to determine whether the formulation of the problem using the "flux" or "com- 
plex” method has any computational advantages in the sonic region over the property- 
derivative method. The perfect gas model also permits a direct comparison of the 
flow-field results with other numerical procedures which use the same model; hence, 
differences in the results cannot be attributed to differences in real gas models. 

SYMBOLS 

Aj coefficients in governing differential equations, defined in appendix A 

Bj coefficients in governing differential equations, defined in appendix A 

b 0 ,c 0 defined by equations (22) 

C coefficient in governing differential equations, defined in appendix A 
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^ref 
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q R 

q Rx’ q Ry 
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coefficients in governing differential equations, defined in appendix A 

functions appearing in governing partial differential equations which are dif- 
ferentiated with respect to y (see eq. (5)) 

static enthalpy 

specific enthalpy divided by gas constant for cold air, °K (see appendix C) 

2 2 

total enthalpy, h + - — * v - 

Ci 

functions appearing in governing partial differential equations which are 
differentiated with respect to x, defined by equations (9) 

nonhomogeneous functions in transformed governing partial differential 
equations, defined by equations (9) 

integers 

Mach number 

&r\ 

pressure 

reference pressure used for appendix C, 1.01325 X 10^ dynes/cm2 
curvature of reference surface, 

R b 

radiation flux vector 

radiation flux components in x and y directions, respectively 

R'T' 

nondimensional gas constant, — — 2L 

WU’oo 2 

radius measured from axis of symmetry of body (see fig. 1) 
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Rjp 
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w ref 
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5 0 


universal gas constant 

gas constant used in appendix C and based on P^ ef , p' re p and W re f, 

T?' — 1 Ot<. — 1 

ref 273.15 

body radius of curvature (see fig. 1) 
radiation flux term, defined in appendix A 
body radius of curvature at x = 0 
temperature, °K 
free-stream velocity 

velocity components in the x and y directions, respectively 
resultant velocity 


mean molecular weight of equilibrium air, g/g-mole (see appendix C) 

mean molecular weight of cold air, 28.96 g/g-mole (see appendix C) 

coordinate along body surface (see fig. 1) 

coordinate normal to body surface (see fig. 1) 

W re f 

compressibility, — — 


W 

axial coordinate (see fig. 1) 
constant 

ratio of specific heats 

shock displacement distance (see fig. 1) 

shock displacement distance at axis of symmetry (x = 0) 
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converged value of shock displacement distance at axis of symmetry 


expression defined in appendix C (eq. (C8)) 
transformed y-coordinate, ^ 
body inclination angle (see fig. 1) 
metric coefficient, 1 + Q6 tj 
expression defined in appendix C (eq. (C16)) 
expression defined in appendix C (eq. (C18)) 
density 

reference density used in appendix C, 1.29313 x 10 "^, g/cm.3 
expression defined in appendix C (eq. (C12)) 
shock-wave angle (see fig. 1) 

Subscripts: 

b refers to body 

j refers to a particular governing differential equation: 

1 shock geometry equation 

2 continuity equation 

3 x-momentum equation 

4 y-momentum equation 

5 energy equation 

s shock-oriented properties (see sketch (a), appendix D) 

o refers to conditions at axis of symmetry on the body 

0 refers to conditions at body surface, 77 = 0 
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1 refers to conditions immediately behind shock wave, tj = 1 

stag refers to stagnation -point conditions 

rj variation with 77 

00 dimensional free-stream conditions 

Primes have been used to denote dimensional quantities. (See eqs. ( 7 ).) Unprimed 
quantities are dimensionless and those with circumflexes p, p, and h are based on 
reference values as shown in appendix C. Double subscripts refer to first, j 
(governing equation number) and second, 77 (location within the shock layer, 0 for 
body and 1 for shock). 

ANALYSIS 


Governing Equations 

The conservation equations for the steady axisymmetric flow of an inviscid radiating 
gas are: 


For continuity: 


For x-momentum: 


For y-momentum: 


For energy: 
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( 1 ) 

(2) 

( 3 ) 

( 4 ) 


These equations have been written in the "divergence” form (ref. 8): 


8Lr aGjfcr 

-^(x,y;ui . • . u k ) + (x,y; Ul . . . u k ) + Kj^x^jU! . . . u k ) = 0 

(j = 1, 2 . . . k) 


( 5 ) 


where x,y are the independent variables, u l • • • % are the unknown functions, and 
Ij,Gj,Kj are the known functions of x,y;u 4 . . . u^. 

The relationship between the shock and the body geometry in the body-oriented 
orthogonal coordinate system shown in figure 1 is 

^ - (1 + Q5)tan(<n - 0 b ) = 0 (6) 


The quantities have been nondimensionalized as follows (primes denote dimensional 
quantities): 
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It is desirable to transform the independent variables x,y to the normalized 
variables x,p where 77 = y/ 6 and 6 = 5(x). When the transformation is made, equa- 
tion (5) becomes 


where 
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The relation K 5 has been established by making the local-tangent slab approximation 
where flow properties for the radiation computations are assumed to be constant in planes 
tangent to the body and therefore 

q Rx “ a^( q Rx) “ ® 
q R = q Ry 

Application of the method of integral relations.- The method of integral relations 
reduces the system of nonlinear partial differential equations to an approximate set of 
nonlinear ordinary differential equations which can be integrated numerically. 

The basic approach is first to divide the region between the body and the shock into 
N strips by using arbitrary values of 77 and then to represent certain combinations of 
the properties by interpolation polynomials in 77 , the coefficients of which are evaluated 
at the strip boundaries. Finally, the system of equations is integrated from 77 = 0 to 
the boundary of each of the strips to give N(j - 1) + 1 independent equations. (It is noted 
that the j = 1 equation relates the shock shape to the body geometry and thus is inde- 
pendent of the number of strips taken.) 

The one-strip approximation is employed in this analysis. For simplicity, the 
functions Ij are assumed to vary linearly in 77 across the shock layers; that is, 

Ij =Ij 0 +(lji -Ijo) 1 ? 0 = 2,3, 4, 5) (10) 

Higher degree interpolation polynomials have been tried in the one-strip approach; they 
complicate the method without a noticeable improvement in the accuracy of the 
computations. 

The governing equation (8) can now be integrated once over the shock layer from 
the body (77 = 0) to the shock (77 = 1) and then differentiated with respect to x to yield 
ordinary first-order differential equations of the form: 

A )f +B j5T + c ^r + E i = ° <i = 2 ’ 3 ' 4 ’ 5 > < n > 

A detailed development of the governing equation (11) and the coefficients are given in 
appendix A. 

Equations (6) and (11) constitute a system of five governing differential equations 
with the five unknown dependent variables 8, u>, I2Q, I3Q, and Igg. (Note that I40 = 0 
for all x when the boundary conditions Vq = 0 and dvg/dx = 0 are applied.) The 
derivatives which are specified by the governing differential equations are computed suc- 
cessively in the following order: 
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(12) 


H = (1 + Q6)tan(w - 6^) 

d„ A 4j^ E 4 

dx B 4 

^20 = A 2 ^ + B 2 Hx + E 2 

dx C 

dI 30 _ ^20 

dx U ° dx 

<**50 _ _ A 5 + B 5 d5T + E 5 

dx C 


(13) 


(14) 

(15) 


(16) 


Note that the exact x-momentum equation (15) is used in lieu of the approximate 
relation. Experience has shown that numerical instabilities can be a serious problem 
in the flow-field calculations. For the one -strip approach with a linear approximation 
(the method being reported), this problem can be eliminated by replacing the approximate 
x-momentum equation at the body with the exact x-momentum equation Equation (15) 
can be used for nonequilibrium and radiating flow-field analyses where the flow along the 
body streamline is not isentropic. 


Initial values .- In order to begin the integration of the governing differential equa- 
tions (eqs. (12) to (16)), it is necessary to obtain initial values of the quantities appearing 
in the right-hand side of these equations. Direct substitution of the initial values at 
x = 0 results in indeterminate (0/0) expressions. Application of L'Hospital's rule yields 
the proper starting values for the nonzero derivatives 


den 

dx 


, 3 + 36 n + 6 


o 2 )|_P 0 " (Pi + Pl v l 2 )_ 


6 0 (3 + 26 0 JpjVj 


(17) 


^20 

dx 


6 0 (3 H- 26 0 ) Pl + (3 + 3 C q + 6 0 2)p 1 v 1 

6 0 (3 + 5 0 ) 


(18) 


For the radiating case, the stagnation enthalpy is also unknown. The expression 
for Hq is 
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(19) 


Hr 


= 6p(3 + 25 Q )p 1 H 1 -g^j g) + (3 + 36 0 + + 38 0 R } 

5 0 (3 + 25 0 ^pj + (3 + 36 0 + 6o^)pi v i 


where Rp is the radiation transport from the slab and is defined in appendix A 
(eq. (A12)). In the present analysis the radiation flux term Rp is set equal to zero. 
(For a radiating shock layer, an iteration scheme can be used in which Rp is initially 
assumed and compared with the computed radiation flux at the stagnation streamline.^ 
Appendix B presents a detailed development of the initial value expressions. 

The equations for the initial values (eqs. (17) to (19)) contain 6 0 the shock dis- 
placement distance at the axis of symmetry (x = 0). Since this term is an unknown in the 
problem, it is necessary to have a basis for establishing the correct 5 0 and hence a 
unique flow-field solution. The logic for the uniqueness is based on the existence in the 
solution of the "sonic singular point." (See ref. 1.) The singular behavior in the solu- 
tion provides a basis for iterating to establish the correct 6 0 . This aspect of the prob- 
lem is discussed in detail in the section on the initial-value sensitivity study. 

Flux and property evaluations.- The basic governing differential equations and the 
coefficients have been developed by considering only the fluid mechanics of the flow field. 
In order to solve these equations, it is necessary to evaluate the thermodynamic proper- 
ties properly. The present analysis is conducted for both perfect gases (y = Constant) or 
equilibrium air. The thermodynamic properties for the equilibrium air model are cor- 
related in appendix C from the equilibrium air data of Hilsenrath and Klein (ref. 9). 

The conditions immediately behind the shock are computed from standard Rankine- 
Hugoniot relations for oblique shocks. The detailed shock equations used in this analysis 
are presented in appendix D for the perfect gas and equilibrium air models. 

Perfect gas analysis: In the perfect gas analysis the stagnation properties on the 
body are computed exactly from the following relations: 


^stag ^1 


J stag 


(y + 1) 2 M| 
4yM^ - 2(y - 1) 

Pstag 


y/(y- 1 ) 


p 




( 20 ) 


In this modified method the properties on the body surface (71 = 0) are functions of 
the fluxes and are computed from the following simple algebraic equations: 
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(21) 


where 


21 


Ur 


20 


-b Q ± l/b Q 2 - 4c c 


b ( 


_2 'X i 30 


° r(2R + l) - l 

(r + 1 ) I 20 2 

y(2R + 1) - 1 

_ x 20 
P 0 u 0 


( 22 ) 


(23) 


P 0 = I 30 - P 0 U 0 2 


(24) 


The sign outside the radical in the Uq expression is positive in the subsonic region and 
negative in the supersonic region. 


The linear variations between body and shock values for Ij^ and also for p^u ?? 2 
are employed. Thus, at all other values of 77 across the shock layer, the properties 
are computed from the following relations: 


I477 
x 2 t 7 


u 


P Ur, 
^ Ti 'I 


_ V 


71 " 1277 

l2?7 

u 

^77 

P 77 = ^-377 - Priori* 


(25) 


Equilibrium air analysis: The thermodynamic properties for the equilibrium air 

model are correlated in appendix C from the Hilsenrath and Klein data for air in chemi- 
cal equilibrium at temperatures up to 15 000° K. (See ref. 9.) The work in appendix C 
is presented in the form of a set of correlation formulas for the equilibrium thermo- 
dynamic properties of air. These formulas are used in the present method since their 
use results in a considerable savings in computational time. The air correlation results 
are of the form: 
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where 


P = P(p,h) 
Z = Z(p,h) 


(26) 


h = H - 


u* + v £ 


The temperature, in degrees Kelvin, is given by 


T , . W refP’ 
ZR'p' 


(27) 


(28) 


In the perfect gas analysis the stagnation properties are computed exactly as func- 
tions of M m and y; this computation is consistent with the formulation of the problem. 
However, in the equilibrium air analysis, which contains the pressure correlation, no 
convenient expressions exist for determining the exact stagnation properties. An 
approach is desired for computing the stagnation properties which is consistent with the 
approximations of the basic numerical approach, is extendable to the coupled radiating 
flow problem, and is not unduly complicated. Two approaches which meet these require- 
ments are: (1) an assumption that the stagnation pressure is equal to the total pressure 

immediately behind the normal shock; that is, 

p stag K Pl + |^l v l 2 < 29) 

which is accurate to within 2 percent of the exact stagnation pressure for hypersonic 
free-stream Mach numbers; and (2) use of the limiting form of the x-momentum equa- 
tion in the set of initial value equations (eqs. (17) to (19)), which gives 


'a2 p 


8x 2 /x=0 


du 0 N 

P 0\dx 


|x=0 


1 f^20 


Pq\ dx / Jx=0 


and an assumption of Newtonian form for , that is, 


\ax 2 y x=0 

“j) x=0 " ^( p stag " p °°) ~ ^ P stag 


(30) 


where /3 is a constant. For a spherical body, modified Newtonian theory indicates that 
(3 = -2.0 and modified theory with centrifugal correction indicates j3 = -2.67, whereas 
the pressure correlation of reference 10 gives (3 = -2.5. The numerical results obtained 
from these two assumptions are discussed in subsequent sections of this report. 
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A Newton-Raphson iteration technique is used to obtain the set of properties on the 
body which are consistent with the Ijq fluxes. The distributions of the properties across 
the shock layer are then computed from the following relations: 


_ *4t7 

77 = 1277 


(31) 


H 


_ 

^1277 


(32) 


where the linear approximations for the Ij^ (j = 2,4,5) terms are employed. (See 
eq. (10).) 

As a result of observing the near-linear resultant velocity distribution across the 
shock layer, which was generated by the inverse solutions, the resultant velocity in the 
present equilibrium air analysis is assumed to vary linearly in 77 ; that is, 

v Rt7 = V R0 + ( V R1 “ Vro) 7 ? ( 33 ) 

The tangential velocity component and the static enthalpy expressions become 

1177 = )jvR ^ (34) 


H t 


V 3 
v Rt? 


(35) 


The pressure distribution is then computed from the quadratic expression 

P77 = Pq + 6Qp o u o 2r ? + ( p i " p o " 6 Qpo u o 2 ) t ? 2 


(36) 


The coefficient of the first-order term in equation (36) is obtained from the exact y 
77 momentum equation upon applying the boundary conditions at 77 — 0 that Vq = 0 
9vq^9x = 0. The remaining state properties are computed from the thermodynamic 
correlations: 

P77 = 


or 

and 


(37) 


T 



(38) 


Initial -Value Sensitivity Study 

The blunt -body flow field is of the mixed flow type which requires the solution to 
elliptic equations in the subsonic region and hyperbolic equations in the supersonic region. 
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The supersonic solution introduces no special problems beyond hyperbolic stability con- 
siderations. However, singularities are present in the governing equations in the tran- 
sonic region which require, for uniqueness, the solution of a two-point boundary problem. 
(See ref. 1.) The subsonic solution, in the one-strip approximation, requires an initial 
value of the shock displacement distance 6 0 (at the axis of symmetry) which will satisfy 
a regularity condition at the natural sonic point on the body, the location of which is not 
known a priori. 

In the property derivative formulation of the problem (refs. 1, 2, 5, and 6), the 
singularity in the transonic region is manifested in the equation for dug/dx which has 
a denominator that vanishes when uq attains the sonic value. The basic procedure is 
to refine 5 0 until both the numerator and denominator of the equation for dug/dx 
simultaneously approach zero within some predetermined accuracy constraint. The 
velocity derivative is then extrapolated into the supersonic region from a point where 
the velocity is within a certain percent of sonic velocity (around 95 percent).. 

Xerikos and Anderson (ref. 6) conclude that 6 0 must be refined beyond eight sig- 
nificant figures in order to maintain downstream accuracy. Calculations (unpublished) 
have been made by Jerry C. South, Jr., of the Langley Research Center by using a 
property-derivative integral method. South’s analysis has shown for the sphere at 

= 5.017, y = 1.4, that for a perturbation in the fifth place in 5 0 , the tangential veloc- 
ity derivative behaves erratically for a small distance downstream of the sonic point. 
However, the integral curve of uq was virtually unchanged from the solution generated 
with a 5 0 which was accurate to about eight significant figures. 

In the flux formulation of the problem (refs. 3, 4, 5, and 11), the boundary condi- 
tions, at the sonic point on the body, of finite velocity Uq and a zero mass flux deriva- 
tive (that is, dl2o/d- x = 0j provide the convergence criterion for the proper 5 0 . 

In the present perfect gas analysis, it is required that dl 2 g/dx and the term within 

the surface velocity expression, \/b 0 2 - 4c 0 (eqs. (21) and (22)) reach zero (within some 
specified accuracy criterion) simultaneously. The initial value of 6 0 is perturbed by 
an automated halving mode on the upper and lower limits of 6o until the sonic-point 
convergence criterion is satisfied. The calculation is then continued into the supersonic 
region by a sign change on the radical of the velocity expression (eq. (21)) from positive 
to negative and by the suppression of any negative values of the term within the radical 
which inevitably arise because of the error in satisfying the exact sonic conditions. 

The saddle-point singular behavior of the subsonic solution is shown in figure 2 for 
the flow of a perfect gas, y = 1.4 over a sphere at = 5.017. The surface velocity 
distributions, calculated by the present method, which uses flux derivatives, and by the 
property-derivative method (South's analysis previously mentioned) are shown. The 
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results are for various perturbations in the respective converged shock displacement 
distances. The 6 0 * is the value of the shock displacement distance which the com- 
puter accepts as being within the specified accuracy of the convergence constraints 
^6 0 * = 0.15297 in the present analysis and 0.15723534 in the property-derivative method 
In the present method the requirements for the unique flow-field solution was for dl2o/dx 

and \[b 0 ^ - 4c 0 to be less than 10" ^ simultaneously. This criterion produced an accu- 
racy of 6 0 * to about four significant figures whereas the S 0 * in the property- 
derivative method was refined to eight significant figures. The comparison in figure 2 
of the behavior of the surface velocity distributions obtained by using the two approaches 
shows that the numerical computations are less sensitive to perturbations in 6 0 in the 
flux formulation than in the property-derivative formulation of the problem. 

The continuation of the solution into the supersonic region is shown in figure 3, 
where the surface velocity distribution is plotted for various deviations in the fourth 
place in the initial shock displacement distance. Values of the surface Mach number Mo 
for the 5 0 - 6 0 * = 0 case are shown in the figure. Although the solutions for the three 
values of 6 0 ~ 6 0 * do agree in the supersonic region, figure 3 illustrates that the initial 
conditions are not sufficiently accurate to produce well-behaved results in the transonic — 
low-supersonic regions of the flow. No notable improvement was observed in the tran- 
sonic solution when the accuracy criterion was refined from 10"^ to 10~5 on dl2o/dx 

and ^b G ^ - 4c 0 since both accuracy constraints produced a 6 0 * to about four-figure 
accuracy. 


The study results indicate that the disadvantage in the present flux approach is that 
the method does not conveniently lend itself to further refinements in 6 0 beyond about 
four figures; hence, a very smooth continuation into the supersonic region is not obtained. 
The Runge-Kutta integration step sizes had to be repeatedly halved in the sonic region 


in order to generate the 10 ^ convergence constraints on \Zb 0 ^ - 4c c 


and 


dl 


20 


dx 


simul- 


taneously and it was noted that values in 6 0 ranging from 0.1527 to 0.1531 would satisfy 
these boundary conditions. Further refinements in the accuracy constraints are undesir- 
able in terms of machine time required by the Runge-Kutta numerical integration proce- 
dure and in terms of round-off errors. 


Thus, the apparent advantages in the flux method over the property-derivative 
method is that the subsonic solution can be generated with less refinement in 6 0 (a 
distinct advantage in performing the time-consuming coupled radiation calculation) and 
that the solution can be continued into the supersonic region without extrapolation of the 
derivatives through the sonic point. However, within the limits of the accuracy on 5 0 *, 
the method does not produce well-behaved results in the transonic region. 
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It should be noted that Belotserkovskii (ref. 5) recommends a hybrid procedure, 
wherein the unknown initial parameters are refined, and the solution advanced by using 
the property derivative equations. When the regularity condition at the sonic point is 
satisfied to the specified accuracy, the flux equations are used to pass through and beyond 
the sonic point. 

Figures 2 and 3 are for the perfect gas model, but about the same techniques can be 
used for equilibrium air to generate similar results. It should be noted that no convenient 
algebraic expression is available for uq for the equilibrium air model and thus the 
analysis requires either an extrapolation of uq into the supersonic region or the selec- 
tion of the proper values of the surface properties which will yield velocities greater than 
sonic just beyond the sonic point. 

RESULTS AND DISCUSSION 

The results of the numerical approach to the blunt-body flow-field solution are 
presented here. The solutions are applicable for the inviscid adiabatic steady flow of a 
perfect gas or equilibrium air over a spherical body. 

Flow-field results are presented for the perfect gas model at a free-stream pres- 
sure of 0.01 atmosphere, a density of 1.562 x 10~® gram per cubic centimeter, and a 
specific heat ratio of 1.4 at free-stream Mach numbers of 5.017, 10.0, and 16.6. The 
equilibrium air cases which are examined are as follows: 

Case I: ulo = 4.572 x 10® cm/sec (15 000 ft/sec); altitude = 45.72 km (150 000 ft) 

Case II: U'^ = 9.144 x 10® cm/sec (30 000 ft/sec); altitude = 45.72 km 
(150 000 ft) 

Case III: U^ = 1.3714 x 10® cm/sec (45 000 ft/sec); altitude = 60.96 km 
(200 000 ft) 

A resume of the free-stream conditions and the results for each of the perfect gas and 
the air cases are given in tables I and II, respectively. 

Flow-field results are compared with the results of the inverse flow-field programs 
of Marrone (ref. 12), Garr and Marrone (ref. 13), Lomax and Inouye (ref. 14), and Inouye 
(ref. 10); with the results of a property-derivative method solution obtained by Jerry C. 
South of the Langley Research Center; and with the experimental data of Sedney and 
Kahl (ref. 15). 


Perfect Gas Results 

Property d istributions _along the body surface.- The body surface pressure distribu- 
tion for the = 5.017 case is shown in figure 4 in both the subsonic and supersonic 
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flow regions. The pressure distributions predicted by "modified" Newtonian theory with 
and without the centrifugal correction are plotted for comparison. The two modified 
Newtonian theories tend to bracket the data obtained by the present method. In the sub- 
sonic region, the present results for the sphere tend to agree more closely with the 
centrifugal correction theory. Also shown in figure 4 is the pressure distribution which 
was correlated from the inverse-flow-field solution. (See refs. 10 and 14.) As evidenced 
by the figure, the pressures agree very well with the inverse correlation except near the 
sonic point. The disagreement in the pressures in the transonic region (x » 0.7 to 0.8) 
is associated with the limited three to four place accuracy of S 0 * and, as previously 
mentioned, this accuracy is not sufficient to produce well-behaved results near the 
singularity. 

The surface velocity and pressure results of the three perfect gas cases are shown 

in figures 5 and 6, respectively. The accuracy criterion of \]h Q ^ - 4c 0 s 10 - ^ and 
dl 20 /dx = 10 _ 3 was required in the = 10.0 and M m = 16.6 solutions. The pres- 
sures in the subsonic region are of acceptable levels; however, in the supersonic region 
for the = 10.0 and = 16.6 cases, the pressures are obviously too low since 
they become negative on the sphere. The surface velocities are correspondingly too high. 
It is further evident in figure 6 that the surface pressures obtained in the present analysis 
are lower in the supersonic region than the experimental values (ref. 15) for the sphere 
at = 5.017 and are lower than both the one -strip property-derivative results 

= 5.017, y = 1.4 and the two-strip results of Belotserkovskii (ref. 5) for Moo = 10.0, 
y = 1.4. 

Shock and bod y shap es.- The stagnation-line shock-displacement distances which 
were obtained from the three perfect gas cases are shown in figure 7 as functions of the 
free-stream Mach number. Also shown in the figure are the results of the two-strip 
property-derivative integral method (ref. 16), the results of the inverse solution (refs. 12 
and 13), the correlation equation given in reference 10, and the experimental data of ref- 
erence 15. The shock-displacement distances computed by the present method are within 
5 percent of the other results presented in this figure. 

Shown in figures 8, 9, and 10 are the shock shapes which were obtained in the pres- 
ent approach for M^ = 5.017, 10.0, and 16.6, respectively. For the M^ = 5.017 case 
(fig. 8), a comparison is shown between experimental data (ref. 15) and analytical results 
of the present method and the inverse method (ref. 13). These data are in good agree- 
ment; however, some deviations begin to appear in the transonic region (x « 0.6). For 
the M^, = 10.0 and 16.6 cases (figs. 9 and 10), a comparison is shown between results 
of the present method, results of the two-strip property-derivative integral method 
(ref. 16), and results of the inverse method (ref. 13). These data also compare favorably 
but indicate some deviations beginning in the transonic region. 
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Property distributions across the sh ock l ayer.- Comparisons of the properties 
across the shock layer which were obtained from the present analysis with the results of 
reference 13 are shown in figures 11, 12, and 13 at various body stations in the subsonic- 
transonic regions. 

Shown in figures 11(a), 11(b), and 11(c) are the pressure, density, and resultant 
velocity distributions, respectively, across the shock layer for M m = 5.017 at x = 0.04, 
0.28, 0.52, and 0.76 radian. The flow is subsonic at the x = 0.52 and lower and is in the 
transonic region at x = 0.76 (Mo = l.l). The property distributions across the shock 
layer agree with the inverse solution to within 2 to 3 percent in the subsonic region. At 
x = 0.76, the present solution gives properties near the body (?y = 0 and 0.2) which deviate 
from those obtained from the inverse solution by approximately 10 percent. This dis- 
agreement is associated with the sonic singularity problem as mentioned previously. 

The property distributions across the shock layer at M to = 10.0 and 16.6 are 
shown in figures 12 and 13, respectively, for various locations in x. As in' the 

= 5.017 case, the pressure, density, and resultant velocity are within 2 to 3 percent 
of the inverse solution values in the subsonic region. In the transonic-supersonic region, 
x = 0.76 and x = 0.73 at = 10.0 and = 16.6, respectively, the body surface 
properties are within about 10 percent of the inverse -solution values. 

Since the inverse solution is restricted to the subsonic -transonic regions of the 
flow field, no comparisons of the properties distributions were made in the supersonic 
regions with this method. 

Comparisons of the density and pressure results which were obtained by the present 
approach in both the subsonic and supersonic regions with the experimental results of 
reference 15 are shown in figures 14 and 15 for the - 5.017 case. There appears 
to be some disagreement in the densities obtained from the experimental data (ref. 15) 
and the present method as shown in figure 14. However, it should be noted that the shape 
and location of the constant density curves are extremely sensitive to the density values, 
particularly in the stagnation region. ^Observe the character and location of the curves 
for p'/p^ - 5.0 and p'/p M = 4.8 where the difference in the density is only 4 percent.^ 
The largest discrepancies in both the densities and the pressures occur near the body in 
the supersonic region where the present approach predicts values too low. 

Belotserkovskii and Chushkin (ref. 17) recommend in the free-stream Mach num- 
ber range of 4 to 6 that a second-order approximation to the fluxes is required for an 
adequate subsonic solution over the blunt body, and above = 10 the first-order 

approximation is sufficient. The present analysis indicates that the first-order approxi- 
mation yields adequate subsonic flow-field results even for the M m = 5.017 case. 
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It is appropriate to note that apparently, the results which one obtains from the 
integral-relations approach depend somewhat on the choice of the quantities which are 
assumed to vary linearly across the shock layer. For instance, Xerikos and Anderson 
(ref. 18) indicate that the use of the unaltered continuity equation, as was employed in the 
present analysis, produces surface pressures which are lower than those obtained from 
a combined continuity-entropy equation, as was employed in both of the Belotserkovskii 
analyses (refs. 5 and 16), and lower than experimental data. This effect was not fully 
investigated in this analysis because the use of the equation for isentropic flow on the 
body surface ( p/pY = Constant) is not convenient in an equilibrium air analysis and is not 
valid for nonadiabatic shock layers (radiating flow fields). However, it was observed in 
the present analysis, that when the exact x-momentum equation was replaced with the 
isentropic expression ( p/p Y = Constant), the results obtained from the two approaches 
were identical. 

Xerikos and Anderson (ref. 18) recommend consideration of the continuity-entropy 
formulation in equilibrium air analyses; however, this procedure requires an additional 
correlation (for entropy) and further it is unlikely that this method can be effectively 
applied to radiating flow-field analyses where the entropy along the body surface is not 
constant. 


Equilibrium Air Results 


As mentioned in the analysis section, there are several ways in which 
stagnation pressures can be generated. The two methods considered here are 


(1) P 


stag 


P 



stag $ \9x 2 /x=0 


Pi + ^Pi v i^ approximation and (2) the second derivative approximation, 


In the subsequent sections on the equilibrium air analysis unless it 


is otherwise specified, the first approximation for P s ^- a g is employed. 

Property distributions along the body.- The body surface pressure distributions for 
the equilibrium air solutions of cases I, II, and III are presented in figure 16. The pres- 
sure distribution which was correlated by Inouye (ref. 10) from the inverse flow-field 
solution of reference 14 is also shown in figure 16. The surface pressures for the three 
cases are within 8 percent of the correlation. 


The surface density distributions for the three cases are shown in figure 17 along 
with the corresponding density distributions of the inverse solution (ref. 14). Of the 
three cases examined, a maximum disagreement of 6 percent occurred in the density 
obtained by the two approaches. 

The surface temperature distributions which were obtained in the present analysis 
are shown in figure 18. The surface pressure results of the inverse solution (ref. 14) 
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were used as inputs to generate surface temperatures in the equilibrium air program of 
reference 19. These temperatures are also shown in figure 18. The temperatures 
generated by the present method are within 2 percent of the temperatures obtained by the 
free -energy minimization approach of reference 19. 

Shock and body shape s.- The initial values of the shock-displacement distances for 
the three cases were found by iterating until three or four significant figures were 
obtained for the converged value of the standoff distance 5 0 *. 

The values of 6 0 * for the three cases are given in table I and are plotted against 
Poo /Pl in figure 19, along with the correlation equation of reference 10. The initial 
shock displacement distances which were generated by the present method are in fair 
agreement, within about 6 percent, with the shock displacement distances obtained by 
the inverse solution for cases II and III, as shown in table II. The case I result for the 
initial shock displacement distance is 5 percent lower than the displacement distance pre- 
dicted by the correlation equation of reference 10 and is 9 percent lower than the value of 
reference 19. 

Shown in figures 20, 21, and 22 are the shock shapes which were obtained in the 
present approach for cases I, II, and III, respectively. The results compare favorably 
with the shock and body-shape results of the inverse solution (ref. 14). 

Property distr ibutions across the shock layer.- The property distributions across 
the shock layer at various body stations in the subsonic region are shown in figures 23, 

24, and 25 for cases I, II, and III, respectively. The pressure, density, and resultant 
velocity distributions of the inverse solution (ref. 14) are also presented for cases II 
and in. A comparison of the results indicates that the pressures and densities obtained 
from the two approaches agree to within 6 percent. The resultant velocities across the 
entire shock layer are in excellent agreement for the two cases. 

Also shown in figures 23, 24, and 25 are the stagnation streamline (x = 0) tempera- 
tures which were produced by the equilibrium air program (ref. 19). The temperature 
results are within 2 percent of the values generated by reference 19. 

Stagnation pr essure assumption.- One of the primary advantages of using the second 
derivative approximation for the stagnation pressure rather than the postshock total- 
pressure approximation is that the former gives a more accurate value of the stagnation- 
point velocity gradient which is required to calculate convective heating rates. Kuby, et 
al. (ref. 7) examined the stagnation-point velocity gradients predicted by the property- 
derivative method of integral relations and concluded that the method consistently pre- 
dicted values which were excessively high. The present method predicts stagnation-point 
velocity gradients which are too high for the postshock total-pressure approximation 
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du 


— ) = 0.69 for case III whereas the second-derivative approach yields values which 

dxjx-0 J 

' 'du' 


are closer to Newtonian and are more accurate I = 0.39 for case III . 

\Wx= 0 


The large d iff erences in stagnation-point velocity gradient however have little effect 
on the other flow-field results. The pressure distributions obtained by two different 
stagnation-pressure approximations are shown along with the inverse -solution pressure 
correlation (ref. 10) in figure 26 for case III. The pressure distribution is slightly 
improved when the postshock stagnation-pressure approximation is eliminated in lieu 
of the second-derivative approximation for /3 = -2.5. The values of S 0 * for the three 
cases increased about 4 percent for the latter approximation and thus improved slightly 
in comparison with the inverse-solution results. About the same percentage changes 
were noted in the other fluid dynamic and thermodynamic properties. 


CONCLUDING REMARKS 


A modified method of integral relations approach for a first-order approximation 
of the fluxes has been used to study the inviscid adiabatic blunt-body flow field. Perfect 
gas and equilibrium air models are considered in the analysis of the flow over spheres. 

The study results indicate that the modified method of integral relations produces 
subsonic solutions in which the shock displacement distance, the shock shape, and the 
thermodynamic and fluid dynamic properties throughout the shock layer are in good 
agreement with the results of the inverse solutions and with experimental data. 

The solution generated near the sonic singularity is somewhat less sensitive to the 
accuracy on the initial shock displacement distance 6 0 in the flux formulation than in 
the standard or property-derivative formulation of the method. The results of the 
perfect-gas analysis demonstrate that in the modified method of integral relations 
approach, a convergent supersonic solution was obtained with about four significant fig- 
ures for the initial value of the shock-displacement distance. However, the present 
approach predicts surface pressures which are low in the supersonic region. Xerikos 
and Anderson's numerical results indicated that the pressures in the supersonic region 
could be improved by employing a combined entropy-continuity equation instead of the 
pure continuity equation as was used in this analysis. In radiating flow-field analyses 
the combined entropy-continuity approach may have merit in the supersonic region; 
however, the additional complications are not warranted in subsonic analyses. 
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The equilibrium air flow-field results indicate that the method can be successfully 
applied to weak radiating flow-field analyses. But, more significantly, the relatively 
uncomplicated method of integral relations is sufficiently accurate and versatile so that 
it continues to offer promise for extension to coupled radiating flow-field analyses. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., July 18, 1969. 
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APPENDIX A 


DEVELOPMENT OF THE GOVERNING DIFFERENTIAL EQUATIONS 


The general governing differential equations which were developed in the analysis 
section (eq. (11)) were written in the form: 


d5 


do; 


— + Bn — — 

dx ■* dx 


+ C 


dl 


JO 

dx 


+ Ej = 0 


(3 = 2, 3, 4, 5) 


These relations were obtained by integrating the governing equation (eq. (8)) over the 
shock layer and then performing the necessary differentiations with respect to x. 

For the one -strip approximation, equation (8) is integrated once from the body 
(77 = 0) to the shock (77 = 1) to yield: 




Kj dT? = 0 


(Al) 


Interchanging the order of integration and differentiation in the first term of equation (Al) 
and integrating the second term by parts gives 

l ( 5 £ ^ - afWj + ( G i Kr )o + 6 j'o k j d ” = 0 (A2) 

Examining equation (A2) it is seen that integrals of the form f Lr dT 7 appear where 

J 0 J 

r = r^ + ^6 cos 0^77 (A3) 

and from equation (10) 

I i = I iO + ( I ji - ho) 71 (3 -2,3,4, 5) 


Note that the r term is treated separately from the linear flux approximation. In 
some analyses (refs. 1, 2, 5, 6, and 18), the Ij and r are "lumped” together as one 
linear approximation across the shock layer ^that is, Ijr = Ijo r b + ( I jl r l “ I jO r b) 1 ^- 

Now defining 

p i * £ ‘i 1 d ” 

substituting equations (A3) and (10) into the preceding equation, and performing the quad- 
rature yields 
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p» _ £*> 

J 6 

Now, defining (see eq. (A2)) 


3 + 


6 cos 6u 


IjO + 3 + 


25 cos 0 b \ 


P 1 




(A4) 


and expanding equation (A4) gives 


s( p i) = I l( 35r bO * iji) + s 2 ™ 6 Obfto + 2 I Ji) 

Differentiating the right-hand side of equation (A5) gives 


(A 5) 


dx 6 


. , 25 cos 0 b , „ 

3 ( I jO + T jl) + fiO + 2 I jl) 


— + 6(3 + 
dx l 


26 cos e b \dlji 


L b 


dx 


. J 0.6 cos 0 b \dIjO 36 a x dr b 62sin 9 b, 0 j ^ d0 bl ( * fi v 
\ r b )~te V j ° W dT " —% "( j ° jl ) d^J (A6) 

The first term in equation (A2) can now be replaced by its equivalent d^Pjjy^dx, and 
becomes: 

^ - Iji(r b + 6 cos + G jlKl (r b + 6 cos 9^ - G j 0 r b + 6 ^ Kj dr? = 0 (A7) 

Substituting the expression for d^Pjydx (eq. (A 6 )) into equation (A7), evaluating the 
remaining terms in these expressions, and combining like terms yields 


T N 26 cos 0 b/T . 

^JO " X jl) + r^ ( X jO " ^l) 


I +63+ r 


26 cos + / + 5 cos e b \dl j0 


dx 


L b 


dx 


q* . dr b 5 2 sin 0 h d0 h / 6 cos 0A 

lf(ho - in) dT — + 2I ») *r + ^ + * -*bT H 1 ' 6Gi0 


+ §£ l 1 Kj dl) = 0 
r b J 0 J 


(A 8 ) 


Equation (A 8 ) can now be written in the desired form: 
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when it is noted that 

£il = ^jldco + ^il^b 
dx 8o) dx 80t, dx 

thus, the coefficients of the governing differential equations become 

26 cos 0 b 


(A9) 


>> 


Aj = 3 


, \ 4 0 COS VU/ \ 

(ho - hi) + — - hi) 


Bj = 613 + 


C = 6 3 + 


26 cos 0b\ 8I jl 

ft 

6 cos $ b \ 


L b 


> (A10) 


Ej = |Sft. 

J r b 


fr T \ dr b J 0 cos cm b osin» b 

?*> + hi) to - + 63 + —%~kk *r - —^10 + 2 V) 


d^ 

dx 


, / 6 COS 6y\ fiK pi 

6(1 + Q5)M + jGji - 6Gjg + — J Kj dp 


J 


G 6 / * 

An examination of the Ei coefficients reveal that terms of the form — \ dr? 

J r b J 0 J 

appear. These terms may be evaluated from the Kj expressions of equation (9): 

k 2 = o 


/drw d0,\ 

k 3 = -P^ - ^ sin e b &r) + Q ( r b + 6 cos V) : 4 

K4 = -P^ b Q + cos e b + 26 Q cos e b p) - Q^r b + 6 cos 6 ^rfjpu“ 


K s ■ I i^iR) 


Replacing p with its equivalent I3 - pu^ gives 

/dr, dfl, 


/ o\/ ar h d£ h \ 

K 3 = -(13 - PU 2 )^ - 6rj sin 0 b —j + Q(r h + 6 cos 0^pu 


K4 = -Ig^r^Q + cos 0 b + 26Q cos 6-^rij + cos 0 b (l + 6Qp)pu2 
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The Kj expressions can be immediately integrated upon assuming a linear variation 
for pu2 across the shock layer. The Kj terms yield 


65 C 1 Tjr _ 35/ t , t \ dr b , S 2 sin 6-^^ ^ ^d 6 h f 

_ j o K 3 d V - - -(I 30 + I 31 )— + r b ?30 + 2l 3l)^ + SQf 3 


26 cos eA 


Ml 


35/ 2 2\ dr b 6sin0 b ( 2 n 2\ d0 b 

¥7 [Pq u o + p i u i ) dT ^ ( p o u o + 2p i u i ) dF 


£ \ - K 4 d, - -3(Q6 + (I 30 + I 3 l) - ^^(l30 ♦ hi) 


6 cos 0 b 


26 Q cos 0 b 


r b d 0 'V" ' r b 


g . C0Se b L 1 i,2 xn1 ,2U. 5 . . . QC0Se b / n , 1 2 


'( p o u o 2 + p i u i 2 ) + 


( p o u o 2 + 2p l u l 2 ) 


Thus 


36/ N dr b S L 26 cos p b\ aI 21 dP b 6 sin % 0 x d0 b 

e 2 - rTp20 + hi) dF + 6 3 + rT IpT dF " rT ?20 + 2I 2l) dT 


b 


/ 5 cos 0A 

+ 6(1 + 6Q)(l + Fp) G 2 1 


Eo = 6 3 + 


26 cos 0 b \ 9I 31 d0 b 


6 cos 0A 




■ b / 30 b dx 


— - + 6(1 + 6Q) 1 + ^I 41 + 6Q 3 + 


26 cos 0A 


~ M1 


36/ 2 2\ dr b 6 2 sin 0 b / 2 2 \d0 b 

rr( p o u o + p i u i ) d^ FT ( p o u o 2 + 2p l u l j dF 


36 T dr b K / 0 26 cos 0 b ^ 3I 41 d0 b 26 sin 0 b d0 b 

E 4 = ^ r 41 dF + % 3 + ^ dF 141 dF 


/ 6 cos 6\\ / v v^v_/o v t-v \ . * 

+ 6(1 + 6Q)fl + — — P j G 41 - 6G 40 - 3 f 5Q + — D j (l 3Q + I 31 ) 


6 cos 0 b \ 


(Equations continued on next page) 
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e 5 


where 


26 2 Q 


cos 0 b , OT \ 36 cos 0 b / „ 

^ (130 + 2I 31 ) + — (p 0 u 0 2 + p lUl 2) 


(All) 


+£ w k ( poV+2piUi2 ) 

36 h . t \ dr b e2sin e b/ T . OT \ de b , J 9 , 26 cos 0 b\ 9I 51 de b 

= r^50 - ISI)^ - — r^— p50 + 2I 5l) dT + 8 1 3 + r fa ' dF 


+ 61 + 


6 cos 0 b 


(1 + 6Q)Ggj + 6Rp 


Rp = (1 + 5Q) 1 + 


6 cos 0 b \ 


iR,! 


- q 


R,0 


(A 12) 
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DEVELOPMENT OF THE INITIAL VALUES 


In order to begin the integration of the governing differential equations (12) to (16), 
it is necessary to obtain initial values for all the quantities appearing in the right-hand 
side of these equations. Direct substitution of these initial values at x = 0 results in 
indeterminate (0/0) expressions. At x = 0, the symmetry conditions are 


_^_9£_8v_d6_ q 
dx dx dx dx 

The nonzero derivatives are 

dco 

dx 

dup 

dx 

In the limit as x approaches zero 


cos e h dr 1 

lim — = lim -t — = lim Q = 1 

x— 0 r b x— 0 ^ x-0 > 


d0 b 

lim — ^ 
x— 0 dx 


-1 


J 


(Bl) 


(B2) 


The limiting form of the coefficients in equation (21) becomes 
A 2 = A 4 = A 5 = B 3 = Ej = E 3 = 0 

A3 = (3 + 26 0 ^Pq - p/j 
B 2 = 6 0 (3 + 25 0 ) Pl ^l 

B 4 = 6 o ( 3 + 26 o)Pi v i -^7 
b 5 = 6 o( 3 + 26 o)Pi h 1 -9^ 
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E 2 - 5 o( 3 + 5 o)Po + 6 o( 3 + 25o)Pi 

E 4 = 6 0 (3 + 26 0 )p 1 v 1 ^ + 2(3 + 36 0 + 5 0 2 )(p Q “ Pi) + ( 6 + 96 o + 46 0 2 )pi v i 2 
E 5 = 6 0 (3 + S o )p 0 H 0 ^ + 5 0 (3 + 25 0 )p 1 H 1 ^ + (e + 95 0 + 46^^11! + 66 0 R F 


where 

du l _ _^1 dw + = _^1 dw 

dx duo dx 30^ dx Ou dx 


From the y-momentum equation (eq. (13)) 


lim 
x— 0 


dco 

dx 




Substituting these expressions for E4 and B4 yields 

do, = ( 3 + 36 o + 6 o 2 )(p 0 - Pj - Pl v l 2 ) 

dx , , 3ui 

6 0 ( 3 + 25 0 )p 1 v 1 — 

In a similar manner, the continuity equation (14) yields 


dI 20 _ _ 6 o( 3 + 26 o)P 1 “^7HF + ( 3 + 35 o + 6 o 2 )pi v i 
dx 6o(3 + f>o) 


and, the energy equation (16) gives 


H 


0 - 



+ 26 0 )pi h i ^ + (3 + 36 0 + 5 0 2 ) Pl v lHl + 36 Q R F 
6 0 (3 + 26 0 ) Pl ^ ~ + (3 + 35 0 + 6 0 2)p lVl 
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A CORRELATION OF EQUILIBRIUM AIR PROPERTIES TO 15 000° K 

By G. Louis Smith and L. Bernard Garrett 
Langley Research Center 

In order to compute the thermodynamic properties of high temperature equilibrium 
air (or any other mixture of reacting gases) from basic principles, it is necessary to 
specify the pressure and temperature. If these two variables are known, the composition, 
density, enthalpy, etc., may be determined directly. Frequently, in flow-field studies the 
enthalpy and density are given, and it is required to determine pressure and temperature. 
In order to accomplish this calculation rigorously, it is necessary to resort to a double 
iteration whereby through some numerical process one seeks the pressure and tempera- 
ture which corresponds to the given enthalpy and density. This procedure is rather 
lengthy and in flow-field studies, where differential equations must be integrated by using 
the results of the computations, the thermodynamics computations may take an excessive 
amount of computer time. 

This problem can be alleviated by use of correlation formulas expressing pressure 
and temperature in terms of enthalpy and density. These formulas are developed from 
the data of Hilsenrath and Klein (ref. 9) which is used as a standard for thermodynamic 
data. Correlation formulas have been written expressing the pressure, compressibility, 
and temperature of equilibrium air as a function of density and enthalpy, for a density 
ratio P'/p ' re f range of 10"^ to 10 and temperatures up to 15 000° K. The pressure is 
correlated to about 5-percent accuracy, the compressibility about 2 percent to 5 percent, 
and the temperature to about 10 percent, although through most of the range the accuracy 
is better than these values. 

The approach taken is to express the equation of state of the mixture in the form: 

P = P(p,h) (Cl) 

The quantities with circumflexes p, p, and h which are used in this appendix are 
defined as follows: 
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h = 


h'W ref 

R' 


so that for conversion to the nondimensional values used in the text 


P 



P - 
h = 


^ref 


hR' 


WrefU'c 2 


For many purposes, the relation (Cl) is all that is required. If temperature and/or mean 
molecular weight are required, the following relations are used. The compressibility Z 
defined by 

w re f 

Z = (C2) 

W 


is correlated in the form 

Z = Z(p,h) 


(C3) 


The temperature T' (°K) is then given by 


T' = 


P 

ZpR;. ef 


(C4) 


For the correlation for p(p,h), the enthalpy range 0 < h < 600 000° K was divided into 
five regions as follows: 


Region I: 0 < h g 5800° K, or approximately 0 < T' < 1500° K. Within this region 
air is calorically perfect for practical purposes, and 


p = (o. 97513 X 10" 3 )hp 


(C5) 


Region II: 5800° K < h £ 10 500° K; or approximately 1500° K < T' < 2500° K. 

Within this region air is a perfect gas, but vibrational modes are excited and energy 
invested in the formation of nitric oxide is significant. Within this region 

p = (0.345 X 10- 2 )h°- 854 p (C6) 


For the remaining three regions, the phenomena are much more involved and the explana- 
tion for the behavior is not obvious. The formulas which were fitted to the data are 
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Region IE: 10 500° K < h < 35 500° K 


logjO P = ^*955 + log^Q p + (0.1545 + 0.0131 log^Q 

+ 0.016^(2.75 - ?)log 10 p - 0.005^(4 - £)(l + log 1Q p) (C7) 

where 

? = 5 log 1Q h - 20 (C8) 

Region IV: 35 500° K ^ h < 178 000° K 

log 1Q p = 1.565 + 1.036 log 1Q p + 0.668(log 10 h - 4.8) + 1.1675(log 10 h - 4.8) 3 (C9) 

Region V: 178 000° K ^ h < 600 000° K 

logjQ P = -3.015 + 1.036 log^Q P + 0.95 log^Q h (CIO) 

These correlations hold for the density range 

10 -4 Ipl 10 

A comparison of the correlation formulas (C5) to (CIO) with data from reference 9 is 
shown in figure 27. The correlation is accurate to about 5 percent or better. 

The partial derivatives of the pressure with respect to the density and enthalpy are 
as follows: 

Region I: 

= (o.97513 x 10" 3 )p 
9h v 

-5= (o.97513 x 10" 3 )h 

dp v 

Region II: 

= (2.946 x 10‘ 3 )ph“ 0 - 146 
8h 

^ = (3.45 X 10 _3 )h 0 - 854 

dp 


Region III: 


4 = ^fo.1545 + [0.0131 + 0.016(2.75 - 2o]log 10 p- 0.005(4 - 2£)(l + log 1Q pj) 

9h h C x 0 
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jP 

dp 


§fl + c[0.0131 + 0.016(2.75 - ?) - 0.005(4 - 



where 

? = 5 log 1Q h - 20 


Region IV: 


Region V: 


_ P 
8h h 


0.668 + 3.5025(log 10 


h 



= 1.036 £ 

dp P 


-^£ = 0.95 £ 
3h h 


-5£ = 1.036 £ 

dp p 

For the correlation of Z(p,h), it is first noted that in regions I and II air behaves 
as a perfect gas so that 

Z=1.0 (h < 10 500° K) (Cll) 


(or approximately T < 2500° K). For higher values of h it was found that curves of Z 
against logjg h for constant p very nearly coincide when translated horizontally an 
amount dependent on p. This relationship is shown in figure 28, where the abscissa is 


i// = log^Q h - 0.044 logjg p - 0.004^1ogjg pj^ - 3.952 (C12) 

A single curve may thus be used to describe Z. This curve is broken into four regions 
in order to fit low-order polynomials to it. 


Region A : 

For 

h S 10 500 

Z = 1.0 

(C13) 

Region B: 

For 

» p ^ 0.55 

Z = 1.0 + 0.53 4? 

(C14) 

Region C: 

For 

0.55 < 4/ < 

1.3 




Z = 2.0 - 

1.78j/ + 0.21^2 + 1.09f 3 - 0.446 

(Cl 5) 
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where 


4*- 

1 

CO 

yH 

11 

(C16) 

Region D: For 1.3 ^ i^< 2.0 


Z = 3.831 - 5.019£ + 3.41£ 2 + 0.24£ 3 

(C17) 

where 


£ = 1-9 - i/' 

(C18) 


The region if/ > 2.0 is outside the range of reference 9, and was not considered 
in the correlation. The boundaries of these regions are shown in figure 29. The dotted 
line corresponds to = 0. 

The results of the correlation formulas (C13) to (C17) are shown in figure 30, along 
with data from reference 9. The comparison is favorable, the error of the correlation 
being less than 2 percent except in the region 150 000 < h < 250 000° K where a 
5-percent error occurs. 

The temperature is computed from equation (C4) by using the p and Z correla- 
tions. The error of the temperature correlation is about 10 percent in some places, but 
is usually less than this amount. (See fig. 31.) 
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RANKINE-HUGONIOT CONDITIONS 


The conditions immediately behind the shock are computed from the following 
Rankine-Hugoniot relations. In the shock-oriented coordinate system (see sketch (a)), 
these expressions in nondimensional form are: 


Shock 



Sketch (a).- Shock and body-oriented properties. 


For continuity: 

For normal momentum: 


V oo = Pl V S 


P, 


P'ooU'j 


+ V, 


2 _ 


Pi 


777t"~2 + Pi V 


Pio U 


(Dla) 


(Dlb) 


For tangential momentum: 


U oo = u s = COS W 


(Die) 
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For energy: 


Hoc — hco + ^ ~ bj + 


u s 2 + v s 2 


(Did) 


In the body-oriented coordinate system, 


Uj = u s cos 


(w - ^b) + v s sin ( w - 6 b) 
= u s sin(o> - 0^ - v s cos^co - 0^ 


(D2) 


where v^ is defined as positive in the y- or 77 -direction. 

The partial derivatives of the velocity components with respect to the shock angle 
are required in the analysis. They are: 


9ui du s 
3 a) dcD 


cos 


(CD- 


9v l du s 


dvg 

dw 


dv 


sin(o) - 0^ - u s sin(co - 0^) - v s cos^cd - 0^ 


> (D3) 


3o) = "duT sin ( w " eb ) + d^7 COS ( W " 0b ) + Us C ° S ( W ‘ M ' Vs Sin ( W " 0 b)J 


Perfect Gas Analysis 

The perfect gas is characterized by the fact that the ratio of the specific heats y 
is a constant. The conditions behind the oblique shock are computed explicitly from the 
following relations: 


Ug = COS Ct) 


v s = 


1 + y - - g -- sin 2 a) 
Y M 2 , sin ci> 

Cl 


Pi = 


1 ‘ puv'JL y 


+ -KC. 

y + 1 \ 


sin 2 u) - l) 


Pi 


sin 2 cd 

1 + — ^ - M^sin 2 cD 


(Equations continued on next page) 


36 


APPENDIX D 



4 y 

y + 1 



sin w cos 




2 

(y + l)M ao sin o> cos 




(D4) 


Equilibrium Air Analysis 

The expressions in appendix C for equilibrium air thermodynamic properties and 
the Rankine -Hugoniot relations (eq. (Dl)) are sufficient to compute the properties imme- 
diately behind the shock. A Newton-Raphson iteration technique is employed to obtain 
a convergent solution. The derivatives of the properties are computed from the deriva- 
tives of the governing Rankine -Hugoniot equations and the equilibrium air-pressure 
correlations. 


d P 1 
dw 



siA/j 1 %>i\ 

Pi 2 V " p i 8h iJ 



9Pl 

9pi 


(D5) 


dhj 

d<u 


sin 


w cos w 




sin 2 u> d Pl 

Pj 3 dw 


(D6) 


dP-L 

dco 


2 sin cj cos w 


1 - 1 - 
p l 


+ 


sin 2 qj dp l 
p 2 dw 


(D7) 


The expressions for 9p^9hj and Sp^Sp^ are obtained from the equilibrium air- 
pressure correlations and are given in appendix C. 
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TABLE I.- RESUME OF PERFECT GAS SOLUTIONS 
(All quantities are nondimensional; y = 1.4] 


M» 



c 

Pstag 

^stag 

Present method 

(a) 

Refs. 12 and 13 

5.017 

0.9229 

5.448 

0.15297 

0.157235 

0.1455 

10.0 

.9228 

6.169 

.12902 

.135240 

.121 

16.6 

.9209 

6.331 

.12461 

.130691 

.120 


a Unpublished data by Jerry South at Langley Research Center. 


Case 


TABLE II.- RESUME OF EQUILIBRIUM AIR SOLUTIONS 
[All quantities are nondimensional except as notedj 




Free-stream conditions 


Case 

Altitude , 
km 

p' 

■^OO) 

atm 

P'oo ’ 

gm/cm 3 

ULo, 

cm/sec 

I 

45.72 

1.445 X 10 -3 

1.837 XlO- 6 

4.572 X 10 5 

II 

45.72 

1.445 

1.837 

9.144 

III 

60.96 

.2226 

.3152 

13.72 


stag 


I 

II 

m 


0.9577 

.9690 

.9699 


Present method 
T 


Stagnation conditions 


Inverse solution 
(ref. 14) 


stag 


11.28 

16.09 

16.86 


stag’ 

OK 

5 113 
8 964 
13 468 


°o 


0.06568 

.04793 

.04458 


stag 


0.9591 

.9701 

.9708 


p stag 

11.34 

16.29 

17.20 


0.0721 

.0495 

.0475 


Free energy minimization 
(ref. 19) 


Pstag 

0.9659 

.9773 

.9775 


p stag 

11.39 

16.50 

17.29 


T* 

stag’ 

°K 

5 024 
8 826 
13 586 
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Figure 4.- Surface pressure distribution. Perfect gas. 
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Figure 6.- Surface pressure distributions. Perfect gas. 
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Figure 9.- Shock and body shapes. Perfect gas ; = 10.0. 
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Figure 10.- Shock and body shapes. Perfect gas; (VU = 16.6. 
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(a) Pressure. 
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(b) Density. 

Property distributions across the shock layer. Perfect gas ; M, 
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Figure 13.- Property distributions across the shock layer. Perfect gas ; = 








□ 


o 


□ 


Present method : 
O Case 1 
A Case H 
□ Case HI 


Inverse solution (ref.10) 



O A 




Figure 17.- Surface density 
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Figure 20,- Shock and body shapes. Equilibrium air ; case I. 
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Figure 21.- Shock and body shapes. Equilibrium air ; case II. 




Figure 22.- Shock and body shapes. Equilibrium air ; case III. 








(b) Density. 

Figure 24.- Property distributions across the shock layer. Equilibrium air ; case II 
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(c) Temperature, °K. 
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Pstag = 




Inverse solution correlation (ref. 10) 
P/P stag = L0 ” ! ' 25 sin 2 X + 0.284 sin 4 x 
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Figure 26.- Comparison of pressure distribution for different stagnation pressure approximations. Equilibrium air; case III 
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Figure 27 



(a) p = 10. 

re p as a function of h for constant density p. 
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Figure 30.- Concluded. 
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